Indeterminacy, Memory, and Motion in a Simple Granular Packing 
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We apply two theoretical and two numerical methods to the problem of a disk placed in a groove 
and subjected to gravity and a torque. Methods assuming rigid particles are indeterminate - cer- 
tain combinations of forces cannot be calculated, but only constrained by inequalities. In methods 
assuming deformable particles, these combinations of forces are determined by the history of the 
packing. Thus indeterminacy in rigid particles becomes memory in deformable ones. Furthermore, 
the torque needed to rotate the particle was calculated. Two different paths to motion were iden- 
tified. In the first, contact forces change slowly, and the indeterminacy decreases continuously to 
zero, and vanishes precisely at the onset of motion, and the torque needed to rotate the disk is 
independent of method and packing history. In the second way, this torque depends on method and 
on the history of the packing, and the forces jump discontinuously at the onset of motion. 

PACS numbers: 45.70.-n 



I. INTRODUCTION 

A. Motivation 

Static granular packings present many difficulties to 
theorists seeking to understand them. One complication 
is indeterminacy. If the grains are assumed to be per- 
fectly rigid but with friction, one cannot find a unique 
solution for the forces between them; there are usually 
many possible solutions Q, 0, 0, 0, H[ . If one calculates 
forces from particle displacements, indeterminacy disap- 
pears, but in many examples, these displacements are 
tiny fractions of a grain diameter. It is astonishing that 
such tiny distances would change the fundamental na- 
ture of the problem. Furthermore, it is not clear how 
the unique state, determined by particle displacements, 
is related to the set of possible forces arising when the 
particles are rigid. These questions are made more press- 
ing by the widespread use of simulation methods that as- 
sume the particles are rigid, and hence somehow choose 
one value for the contact forces out of many possible so- 
lutions. 

Next, there is the question of memory. The method of 
construction of a packing and its history change its be- 
havior. For example, it has been shown that the method 
of construction changes the distribution of forces under a 
sand pile |fj, |8| • In this case, the memory of the pack- 
ing seems to be stored in the orientation of the contacts 




Finally, there is the onset of motion: what forces are 
needed to permanently deform a static packing, and how 
does this motion originate? These questions are the prin- 
cipal concern of soil mechanics. Motion often is very lo- 
calized, i.e., in shear bands Ilfjj. Another phenomena 
related to the onset of motion is the linear accumula- 
tion of deformation in granular assemblies subjected to 
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FIG. 1: A disk supported by two walls through contacts a 
and p. The angle <j> suffices to characterize the geometry. 

periodically varying forces 0, 0] . 

In this paper, we consider a very simple system 
where the relationships between these problems is clearly 
demonstrated. We hope that it will provide clues on how 
to understand much more complicate packings. 



B. Synopsis 

The granular "packing" shown in Fig.^is perhaps the 
simplest possible granular packing. Nevertheless, it ex- 
hibits many subtleties, and in this paper we use it as an 
example to show the relationship between indeterminacy, 
memory, and motion. This configuration is similar to the 
rod wedged between two converging walls discussed in 
Ref. la. It is even simpler than the one considered in 
Ref . [X3| , where a disk is placed in an asymmetric groove 
and subjected to gravity. Here, we consider only a sym- 
metric groove, but the forces we impose on the disk are 
more complicated. In addition to a gravitational acceler- 
ation g, we apply a torque r > 0. The contacts between 
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the disk and the wall are assumed to be non-cohesive and 
frictional, with Coulomb friction ratio /x. We investigate 
two questions: How do the contact forces at a and (3 
change in response to the imposed forces? At what value 
of the torque does the particle begin to rotate? 

We first examine some general considerations: the 
equations of static equilibrium, and the status of the con- 
tacts. Then we attempt to deduce the contact forces from 
these considerations. We show that for certain geome- 
tries and torques, one cannot decide whether the disk ro- 
tates or not. We th en p resent a series of contact dynam- 
ics (hereafter CD) [lj] and molecular dynamics (here- 
after MD) simulations, and point out the difference 
between them. Then we formulate a second method of 
calculating the contact forces, assuming the arise from 
small deformations of the particles. This second method 
explains all the features of the MD simulations, and sheds 
light on the questions mentioned above. 



It is convenient to use the following orthogonal basis 
of K 4 : 





(5) 



If F is written in this basis: 
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(6) 



then one can verify that 



II. RIGID PARTICLES 



cF = 



a x r cos < 



agr sin < 



= -fe 



(7) 



A. Static Equilibrium 

The equations of static equilibrium for the disk in 
Fig. n are 

R a sin 4> + T a cos — Rp sin <f> + Tp cos 4> = 0, 
R a cos 4> — T a sin 4> + Rp cos 4> + Tp sin </> = mg, 

rT a +rTp = -r. (1) 



This equation enables one to easily construct solutions to 
Eq. |J2j, for one has immediately a x = 0, a y = mg, and 
ag = — t / (r sin 4>). However, ao cannot be determined in 
this way because Fo is a null eigenvector of c. This is an 
expression of the indeterminacy in the problem. 



B. Contact Status 



Here, R a indicates the normal force at contact a, with 
R a > for repulsion. T a is the tangential force, with 
T a > when this force exerts a positive torque on the 
disk. Rp and Tp arc defined in the same way. 
These equations can be written in matrix form: 



Of course, one does not have absolute freedom to 
choose ao, because the following conditions must be sat- 
isfied at the contacts: 

R > 0, fiR> \T\, (8) 



cF + f cxt = 0. (2) 
We call c the contact matrix: 




Note that the dimensions of c require that it have at 
least one null eigenvalue. The contact forces F and the 
external forces f cx t are 




where R a (Rp) are the normal components forces exerted 
at contact a (/?), and T a and Tp are the tangential com- 
ponents. 



where the constant [i is the Coulomb friction ratio. The 
first inequality excludes cohesive forces, and the second 
requires that the tangential forces not exceed a certain 
threshold. 

At this point, it is useful to introduce the concept of 
contact status. If we have equality in the second condi- 
tion in Eq. ©, the contact is said to be sliding, otherwise 
it is non-sliding. The tangential relative motion vt must 
be consistent with the contact status. If the contact is 
non-sliding, no tangential motion is allowed, but if the 
contact is sliding, then the tangential force must oppose 
the motion. Note that v t — is allowed under all cir- 
cumstances, so that contacts in static equilibrium can be 
sliding. 

To see how this limits the possible values of ao, one 
can first use the definition of F in Eq. an d the basis 
in Eq. (J5J to compute the forces in terms of fj,, <j>, and 
the coefficients in the expansion of Eq. ©. Then one 
can use Eq. Q to eliminate all the coefficients except 
for ao- Finally, one forms the inequalities fJ.R a > — T a , 
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\iR a > T a , fiRp > —Tp, and [xRp > Tp. These yield 

((j, - tan<£) > 0, (9a) 
O + tan^) > 0, (9b) 
((j, + tan 4>) > 0, (9c) 
(ju-tan<£) > 0, (9d) 



Some of these conditions are redundant. For example, 
Eq. (|9"c|l implies Eq. (|9b|) . To see this, note that due to 
our choice r > 0, 



2r 



(fx + tan0) > 0. 



(10) 



Adding this to Eq. Jgcj! gives Eq. JSg. Thus Eq. (JM 
never needs to be considered; it will be satisfied as long 
as Eq. (|9"c|) is. Similar reasoning can be used to show 
that Eq. I|9a|) is redundant when [i — tan <p > 0, and 
that Eq. I|9d(l is redundant when /i — tan0 < 0. When 
/i = tan0, these conditions are equivalent. 

The conditions in Eqs. © put limits on the value that 
Oo can attain £|. Furthermore, when a contact is sliding, 
we have equality in one of the above cases, and ao is 
determined: indeterminacy disappears. 



C. Application 

We now attempt to deduce the contact forces and 
particle motion, using only the considerations presented 
above. No relation between particle deformation and 
contact forces is assumed. Note that the biggest diffi- 
culty will be determining the coefficient ao; the other 
three coefficients in the expansion in Eq. © can be read 
off from Eq. Q. The unknown value of ao is an ex- 
pression of the indeterminacy of the problem. We must 
consider three separate cases, depending on the relation 
of jjL to the slope of the sides of the groove, tan0. 



1. Shallow slope: tan0 < /i 

Let us first restrict ourselves to tan0 < [i and fi < 1. 
Under these assumptions, we have /i — tan</> > 0, so the 
ao is constrained by the conditions Eqs. (|9"cl) and l|9d|l . 
Furthermore, /itan</> — 1 < 0, so that Eq. l|9"cl) sets a 
lower, and Eq. I|9d[) an upper bound on ao: 



< a < a„ 



(11) 



where 



mg 



/i — tan (/) 
1 + [i tan ( 



mg 



r sin ( 



[i + tan (f; 
1 — fi tan i 



(12) 



The quantities in the curved brackets are always positive, 
while the quantity in the square brackets is positive for 
t = 0, and decreases as r increases. It vanishes when 



Ti = mgr sin <p, 



(13) 



and becomes negative when r > Ti. When it is negative, 
no solution for ao exists, for a m j n > a max . Thus there are 
no solutions of the equations of static equilibrium, when 
t > Ti that also satisfy Eq. (JSJ, and the disk must move. 
On the other hand, when r < t±, an infinite number of 
such solutions exist, one for each value of ao satisfying 
Eq. (|llfl . Finally, at r = t\ , we have a min = a = a max 
0. hence there is a unique solution. 

We have shown that the disk must rotate when r > t\ , 
and that static solutions exist for r < T\, but we have 
not yet shown that the disk will not move when r < t\. 
After all, for some choice of ao, it might be possible to 
put the disk in motion. But this can be shown to be 
impossible by considering what happens at r = t\. In 
this situation, we have ao = 0, so the contact forces are 



R, 



a = mg cos ( 



-mg sm< 



R 







2> = 0. 



Thus the disk is entirely supported by contact a. R a 
and T a cancel the gravitational force, and t balances the 
torque rT a = — t\ exerted on the disk by the contact 
force. A small increase in r and T a will cause the disk to 
roll upwards to the left, out of the groove. On the other 
hand, if r < t±, then the torque rT a — — n will not be 
balanced, and the disk will try to roll down the slope. 
Therefore, if r < t±, the disk cannot rotate. 

Another possibility is that the disk rotate in place. But 
this is impossible, when tan^ < /z, because this requires 
T a = —^R a and Tp — —fiRp, i.e., contacts a and f3 must 
be sliding. This occurs when we have equality in both 
Eqs. H9a|l and l|^cjl . But this cannot happen because when 
equality holds in Eq. Ij9a(l . then Eq. I|9dj) is violated. 

Therefore, the value of ao has no effect on the motion 
of the disk If t < t\ 1 the disk remains in the groove, 
although the forces cannot be uniquely determined. At 
r = ri, there is a unique solution for the forces: contact 
(3 opens and the particle's weight is supported by contact 
a. For t > T\, the disk rolls out of the groove. 

Note that there is a relationship between the onset 
of motion and the disappearance of indeterminacy: we 
have motion for t > t\ and indeterminacy for r < t\. 
Indeed, we could have anticipated this when we wrote 
down Eqs. Jj|J). For a disk to move, at least one contact 
must be sliding (or open), leading to an equality in at 
least one of Eqs. 10 • Then this equality can be used to 
determine ao, eliminating the indeterminacy. 
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2. Intermediate slope: fj, < tan^ < l/fi 

Now let us consider intermediate angles [i < tan0 < 
At these values of (f>, the relevant conditions from 
Eqs. @ are Eqs. lj§a|l and (jHcf. Eq. (J5cJ sets an upper, 
and Eq. (j§a|) a lower, bound on ao- Thus we have Eq. ltTT|) 
again, except now 



mg + 



r sm < 

T 



rsm< 



tan ( 



1 + fi tan (f. 
tan (f> + [i 
1 — /i tan ^ 



(15) 



Again, the quantities in the curved brackets is positive for 
(i < tan</> < l/fi. The quantities in the square brackets 
are equal and positive for r = 0. But as r increases, a mul 
increases while a max decreases. When 



t = t 2 = mgr- _ 
1 + /i 2 

we have a min = a = a max = a 2 , where 



a 2 = mj 



1 



fi 1 + tan 2 



tan< 



1 + M 2 



tan 4> — 
1 + /i tan <^ 



(16) 



(17) 



and there is a unique solution for the contact forces. 
When t > t 2 , no static solution satisfying the contact 
conditions exists, as a mm > a max . 

Furthermore, we can show that the disk cannot move 
when t < t 2 . The disk cannot roll out of the groove, 
because the forces in Eq. fTty violate \T a \ < fiR a . If 
the disk is to rotate in place, we must have T a = —fiR a 
and Tp = —fiRp, i.e., we must have equality in Eqs. I|9a|) 
and (|9"cj) . But we have just showed that requiring these 
equalities leads to r = r 2 , a m j n = a = a max . 

Therefore, we have exactly the same situation as in 
Sec. Ill C"T1 For t < T2, there is indeterminacy without 
motion, and for r > t 2 there is motion. There is a unique 
solution for the forces precisely at r = t%. Furthermore, 
note that n = t 2 when tan (f> — fi and when tan 4> = l/fi. 



3. Large slopes: tan<jf> > l/fi 

At large angles (tan^ > l//i), the relevant conditions 
are still Eqs. (|9a|) and (|9^|l as in the previous section. But 
now the factor 1— /itan</> becomes negative, transforming 
Eq. l|9*c|l into a Jower bound on a . Eq. I|9a|l remains 
a lower bound on ao, meaning that there is no upper 
bound on ao- It can be made arbitrarily large. Therefore, 
for any r > 0, it is possible to satisfy all conditions in 
Eq. © simply by making ao very large. (Note that all 
the factors multiplying ao in these conditions are positive 
when tan</> > l//x). This means that an infinite torque 
can be put on the disk without causing it to rotate. 

On the other hand, if r = r 2 , one can obtain equality in 
Eqs. I|9a|) and (|9"c|) . meaning contacts a and (3 are sliding. 
If r is increased slightly, then the disk can start to rotate 
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FIG. 2: The f orces applied to the disk. Time is given in units 
of t* — ^/d/g. The gravity attains a maximum at t = tA, and 
the torque begins to increase at t = tB. 



in place. Therefore, both static and moving solutions co- 
exist for r > t 2 , and we cannot decide on the basis of 
Eqs. © whether the disk rotates or not. 

This remarkable situation poses two questions: First, 
how will the CD algorithm handle this case? This algo- 
rithm tries to use the information presented in Sec. ITTIto 
deduce the motion of the particles, but we have shown 
that this is insufficient. Secondly, it would be very easy 
to construct an experiment to measure the torque needed 
to make the disk rotate. One could construct a groove 
with tan</> > l//i, and place a cylinder in the groove, and 
try to turn the cylinder. What determines whether the 
cylinder rotate? These questions are investigated in the 
next sections. 



III. SIMULATIONS 

To verify the conclusions of the previous section, and 
to clarify the situation for tan</> > l/fi, we carried out 
CD and MD simulations. The disk was placed in the 
groove without any forces acting on it. Then gravity was 
turned on, and increased, until it reached a maximum 
of g m ax at t — t A . Then it was decreased, reaching g» 
at t = ig, and thereafter was held constant. A linearly 
increasing torque r was then applied until the particle 
began to rotate. Fig.[2]shows the gravity force and torque 
as functions of time. In the following, we will vary <jf ma x, 
but keep g* fixed. In all cases, [i = 0.5. 

The experiment can be considered as a very simple 
soil mechanics experiment. First, the sample, consist- 
ing of one disk, is "prepared" by pushing the disk into 
the groove, and "loaded" by the torque until it "yields" . 
We are especially interested in knowing if and how the 
preparation affects the yielding torque. 

A typical result of a simulation is shown in Fig. [21 The 
three non-zero coefficients a y , ag, and ao of the expansion 
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FIG. 3: The coefficients of the expansion, Eq. @, as a func- 
tion of time. Here 4> = 20° and \x = 0.5 so that tan0 < (J,. 
The curves stop when the particle begins to rotate. 
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FIG. 4: The coefficient oo as a function of time for various 
values of g ma x/<7*. Here <j) = 20° and = 0.5, so that tan</> < 
jj,. The thin lines are the CD simulations, with g m ax/g = 2 
(lowest curve), 4, 6, and 8 (highest curve). The thick line are 
the corresponding MD simulations, which fall on top of each 
other in this case. The large dot indicates the point where 
the disk begins to rotate. 



in Eq. © are shown. In accordance with Eq. Q, a y and 
ag are identical for the CD and MD solution methods, 
and given by a y = mg, ag = —r/(r sin <f>). The coefficient 
ao, however, differs between the CD and MD solutions. 
In the following, we will plot only ao, as the other coef- 
ficients are always uniquely determined by the imposed 
forces. Three different values of 4> wm be considered, 
corresponding to the three subsections of Sec. Ill CI 



FIG. 5: The undetermined coefficient ao from Eq. @ as a 
function of time. Here, <f> = 50° and \i = 0.5, so that fi < 
tan</!> < l//i. The thin lines are the CD simulations, with 
gmax/ff = 2 (lowest curve) , 4, 6, and 8 (highest curve). The 
thick lines are the MD simulations with the same values of 
Omax/<?- The large dot indicates the point where the disk 
begins to rotate. 

A. Shallow Slopes: tan<jf> < fi. 

In Fig. 01 we plot ao as a function of time for the MD 
and CD simulations. In the MD simulations, ao = 0, 
independent of time and of g max - In the CD simulation, 
ao is proportional to the gravity for t < ts- Note that 
the values of ao are significant - almost twice the weight 
of the particle for <? ma x/3* = 8. However, at t = ts, 
all CD simulations have the same value of ao. As the 
torque is applied, ao moves linearly towards 0. This oc- 
curs because when r < n, the disk cannot move, but ao 
is constrained between the two values given in Eq. 112|l . 
which approach each other as r approaches r±. Indeed, 
the diagonal line segment near 50 < < 60 traced out 
by the CD simulations corresponds to requiring equality 
in Eq. JjjHJ. The conditions Eqs. (|5c|> and (|9dl) act as a 
"funnel" that guides ao toward as the torque increases, 
so that ao = when r = T\. Thus the MD and CD both 
predict that the disk moves at r = T\. However, the two 
methods predict different routes to failure. MD predicts 
that both contacts remain non-sliding until the force at (3 
vanishes, and the disk starts to roll. On the other hand, 
CD predicts that the contact at /3 first becomes sliding, 
and then later starts to roll. 



B. Intermediate Slopes: /i < tan<^> < l//i. 

In Fig.El we show ao for the case of intermediate slope 
/i < tan0 < 1///. The behavior of the CD simulation is 
quite similar to the preceding case. For t < ts, ao is again 
proportional to gravity, but attaining even larger values 
than before. At t = is, all CD simulations are identical, 
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FIG. 6: Coefficients of the expansion Eq. @ as a function of 
time. Here, <j> = 80° and fi = 0.5, so that tan0 > The 
thin lines are the CD simulations, with g ma x/g = 2 (lowest 
curve), 4, 6, and 8 (highest curve). The thick lines are the 
MD simulations with the same values of p m ax/g. The dots 
indicates the points where the disk begins to rotate. The 
diagonal dot-dashed line is obtained by requiring equality in 
Eq. G3. 



and they evolve together during the loading. They again 
encounter the "funnel" arising from the lower and upper 
bounds on ao given in Eq. i|15|) , and move toward clq = a^, 
and the disk starts to rotate when r = t%. 

The MD simulation has a quite different behavior. For 
t < tA, ciq is also proportional to the gravity, although 
lower than the CD value. But for t > tA, ao remains con- 
stant, so that ao is proportional to g max , even at t = ie. 
In this case, therefore, the preparation does make a differ- 
ence. The MD simulations remember the value of g max , 
and this memory is stored in ao- This simple example 
clearly demonstrates a principle stated in the paper's in- 
troduction: the indeterminacy of forces in perfectly rigid 
particles is related to history-dependency in packings of 
deformable particles. The coefficient ao, which is un- 
determined in the perfectly rigid case, here contains the 
memory of the packing, and is simply proportional to the 
greatest downward force that the disk has experienced in 
the past. 

The history dependency of the MD simulation, how- 
ever, does not affect the yielding torque. The different 
MD simulations encounter the funnel at different times, 
but are all guided towards the point r = ra, ao = a2, 
where the disk begins to rotate. Therefore, the memory 
on the MD simulation can only be detected by inspecting 
the contact forces. 



C. Large Slopes: tan0 > 1/fJ,. 

Finally, in Fig. |SJ we show the case where tan <j> > 
1/fx. This case has some important differences from the 



preceding cases. First of all, both the CD as well as 
MD exhibits history-dependence, as the contact forces 
at the end of the sample preparation and the yielding 
torque depend on g max . Secondly, only half of the funnel 
observed at lower angles operates as before. 

One may ask how the CD simulation can exhibit mem- 
ory, because it deduces the contact forces from the prin- 
ciples stated in Sec. [H] but the history of the packing 
does not enter into those considerations. The memory 
of the CD algorithm comes from the way it chooses the 
contact forces. It begins with an initial guess, which it 
then refines through an iterative process until it arrives 
at a satisfactory solution. One usually takes the initial 
guess to be the solution of the previous time step because 
convergence is faster. But in Fig. [S| it is clear that this 
choice of initial guess also serves to make the CD algo- 
rithm history dependent, as the solution chosen depends 
on the initial guess p| . But this memory is not the same 
as in the MD simulations. Nor is it clear why the prepa- 
ration of the sample makes a difference at 4> = 80° but 
not at <j> = 50°. 

But the most important difference with the previous 
examples is that only half of the funnel works as be- 
fore. As was noted earlier, there is no longer an upper 
bound on ao; Eqs. (|9a|l and l|§c|) are both lower bounds. 
From the figure, one can see that the lower bound set 
by Eq. (|9a|l functions as before. When ao reaches this 
lower bound, it then increases linearly as the torque is in- 
creased, until equality is obtained in both Eqs. (|9a|) and 
(|9"c|) . Then the disk begins to rotate. The behavior of 
the lower bound set by Eq. l|§c|) is quite different. When 
equality is obtained in that condition, the value of ao 
jumps discontinuously, and the particle begins to rotate. 
As a result, the yielding torque is determined by when the 
system first satisfies equality in Eq. I|9c|l . and this in turn 
depends on the value of ao set by the preparation phase 
of the experiment. Thus in this case, the memory does 
affect the yielding torque. In Sec. IIV Bl we will explain 
this discontinuity in ao theoretically. But one can already 
understand its physical origin. Equality in Eq. (|9c|) cor- 
responds to contact (3 becoming sliding, meaning that 
tangential motion is now possible there. Since contact a 
is non-sliding, the particle pivots about contact a. As we 
showed earlier, contact (3 cannot open when tan (f> > /a, 
but the disk can apparently move enough to reduce the 
contact forces very quickly, allowing the particle to ro- 
tate. 

In Fig. U| we show the yielding torque as a function of 
<?max- It is initially independent of g max , corresponding 
to the case where the system first meets the lower edge 
of the funnel. Then, CD and MD give different yielding 
torques. The MD values are well predicted by a result 
that will be obtained in Sec. lIVB~3l 

The CD values for the yielding torque are equal to or 
below the MD ones, because ao decreases as the torque 
is increased, while in the MD case, ao remains constant. 
This is plainly visible in Fig. The different CD points 
correspond to different iteration procedures. The CD al- 
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FIG. 7: The torque at which the disk begins to rotate at 
4> = 80° and /i = 0.5, for different values of <7 max . The four 
different sets of points for CD correspond to the different it- 
eration algorithms described in the text. The line shows the 
theoretical prediction given in Sec. IIV B 31 

gorithm searches for a solution by adjusting the contact 
forces one by one. It considers a given contact, and calcu- 
lates the change in its contact forces needed to prevent in- 
terpenetration and minimize sliding. It then adds this to 
change to the its current guess for the forces. After pass- 
ing over all the contacts a certain number of times, the 
changes in the force needed fall below a certain threshold, 
and the solution is accepted. But this procedure can be 
altered by multiplying the calculated change in forces by 
a number A, with A = 1 corresponding to the usual itera- 
tion procedure. Changing A changes the yielding torque, 
as shown above. Furthermore, the memory in the CD al- 
gorithm can be erased at each time step by always using 
F = as the initial guess, instead of the solution of the 
last time step. When this is done, one obtains the points 
labelled "RESET" in Fig. □ and the yielding torque is in- 
dependent of <7ma X - This result emphasizes the important 
role of memory in this experiment. 

IV. DEFORMABLE PARTICLES 

A. Formalism 

We now present a second method of calculating the 
contact forces that is able to explain all features of the 
MD simulations presented in the previous section. The 
cost of this additional information is that more assump- 
tions must be made. First of all, one must specify how 
the forces depend on the deformations, but more signifi- 
cantly, one must specify the past history of the packing. 
This method shows how indeterminacy is replaced by 
memory. More precisely, when tan</> > [i, we show that 
ao in Eq. (JSJ) stores information about the largest down- 
ward force that has been exerted on the disk in the past. 



Indeed, ao is a quantification of the "degree of wedging" 
discussed in Ref. P- 

We model the deformations in a very simple way, as- 
suming that contact forces are generated by springs that 
are stretched by the motion of the disk. Our model is in- 
spired by the MD simulation method |15|. but we apply 
it analytically to our simple problem. A more general 
and realistic calculation has been presented before [l3| . 
but our purpose here is to show how the indeterminacy 
is replaced by history-dependence when forces arise from 
deformations. 

When two bodies touch, a normal and a tangential 
spring be created at the instant of contact. The contact 
forces are simply proportional to the spring lengths: 

R=-kS n , T=-k5 t , (18) 

where S n and St are the normal and tangential spring 
lengths, and k is the spring constant. One normally 
includes damping in Eq. (|18f) . but we will assume that 
the motion is quasi-static, i.e., a sequence of equilibrium 
states. Under this assumption, the particle velocities are 
vanishingly small, and so are the damping forces. How- 
ever, the damping has important consequences that will 
be discussed below. Eq. (jTHjl leads to a vector equation 
for F in terms of the spring lengths D: 

F = -feD, (19) 

where 6 n , a , <5t,a, 5 n ,p, and 5t.p are gathered into D just 
as the contact forces are arranged in F [see Eq. (0J]. 

In general, one could introduce different spring con- 
stants for the normal and tangential springs, and the 
ratio of these spring constants affects the behavior of the 
system, even in the limit of infinitely hard particles. We 
set both spring constants equal, because it simplifies the 
results, and our purpose here is only to show how the in- 
determinacy is removed. But we remark that the ratio of 
the spring constants does not appear when the particles 
are rigid, indicating that the physics associated with this 
parameter has been eliminated. 

The spring lengths obey 

6 n = v n , 5t = {% Vni (20) 

where v n and Vt are the normal and tangential compo- 
nents of the relative velocity. There are two choices for 
5t because fiR > \T\ means that the tangential spring 
has a maximum allowable length: \5t \ < fJ>5 n . If applying 
S t = Vt would lead to a violation of this condition, the 
second choice is taken. Eq. (|20|l relates the time deriva- 
tive of D to the velocity of the disk. If all contacts are 
non-sliding, we have 

3n,a = v x sin 4> + v y cos (f), 

&t.a = v x cos <j) — Vy sin <f> + ruj 1 

S n ,/3 = ~v x sin 4> + Vy sin <j>, 

$t,i3 = v x cos <j) + v y sin 4> + ruj. (21) 
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Gathering the velocities of the disk into a single vector: 




Eqs. (|21[1 becomes 



D 



T 
c 1 V, 



(22) 



(23) 



where c T is the transpose of c in Eq. J3J . The appearance 
of c T is not a coincidence, but a general property valid 
for all granular packings |16| . 

It remains to incorporate the status of the contacts into 
our formalism. This can be be done by inserting a matrix 
S that depends on the contact status into Eq. (12311 : 



D = Sc J 



where 



S a 
S/3 



(24) 



(25) 



Here, S a = 1 if contact a is non-sliding, S a = if it is 
open, and 



S a 



1 
±n 



(26) 



if it is sliding. 

Finally, let us not forget Newton's equation of motion: 



Mv = cF 



(27) 



Here, the matrix M contains the mass m and moment of 
inertia I of the disk on the diagonal: 



M 




(28) 



Eqs. (|T§j) . J23J, and lj27|) form a set of equations that 
can be solved for the motion of the disk. Note that c 
remains constant because all contacts are between the 
disk and the straight walls. The status of the contacts 
can change, however, so it is useful to divide time up into 
segments [to, ti], [ti,t%],..., where the changes of contact 

status occur at the times t\, In the interior of a 

time interval, the matrices c, S, and M are constant, so 
Eqs. (fH?j) . (J2U, and (|2"7|) can be combined into 



Mv = -Qv + f c 



(29) 



where Q = fccSc T . The matrix Q gives the change in 
force that arises from a small displacement of the disk. 
From this equation, the importance of the sign ofv T Qv 
can be seen. If v T Qv > 0, then Q acts like a nonneg- 
ative number in Eq. I|29|l . so that the contacts generate 
forces opposing the velocities v, and hence v will tend to 



oscillate about an equilibrium value depending on f cxt . 
If the damping is chosen properly, these oscillations will 
be damped out very quickly, so v will always be close 
to its equilibrium value. Physically, this means that v is 
slaved to f ox t, as long as f ox t changes much more slowly 
than the damping time scale. This amounts to removing 
the inertia from Eq. I|29|) and obtaining 



Qv 



(30) 



The situation is different if f cxt causes displacements that 
are in the null space of Q, i.e., v T Qv = 0. In this case, 
no forces opposing the displacements are generated, and 
they can grow without bound, leading to permanent re- 
arrangements of the particles. Another possibility is that 
v T Qv < 0. In this case, the contact forces amplify the 
velocities, which then grow exponentially. In this case, 
the packing fails catastrophically, with the contact forces 
changing much more quickly than f GXt . 

Finally, let us note there is a compatibility condition 
between S and v that must be satisfied: v must be such 
that the tangential springs stay stretched to their max- 
imum lengths. If this condition is not fulfilled at a con- 
tact, then that contact becomes non-sliding, and S must 
be modified. 

Let us now show how all this can be used to calcu- 
late the contact forces without indeterminacy. We will 
do so by calculating the coefficients of the expansion in 
Eq. ©. We can construct an equation for these coef- 
ficients by differentiating Eq. (|19|l with respect to time, 
and combining it with Eq. Q24[l . leading to F = Sc v. 
Then we project this equation onto the basis given in 
Eq. (JSJ) by left-multiplying it by a, a 4 x 4 matrix whose 
rows are the basis vectors in Eq. JSJ), but without the 
factor of 1/2. After doing this, we obtain 



(31) 



Eq. (|31|) represents four equations, one corresponding to 
each of a x , a y , ag, and ciq. The first three coefficients 
are known from Eq. J7J), thus the first three equations 
above can be used to find the three components of v. 
Then the last equation can be used to determine do- The 
displacement of the disk from its equilibrium position is 
obtained by integrating v. 




B. Application 



1. Shallow angles: tan<^> < /j, 



When tan <f> < fi, the contacts remain non-sliding until 
the disk begins to rotate. Therefore, we can consider the 
entire experiment to take place during one time interval. 
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When the contacts are non-sliding, 



fcaSc J 



2k 



1 r cos ( 

1 

r sin <; 





(32) 



Using this result in Eq. (|31fl . we have 



he 
do 



2kv x + 2krto cos < 
2kv y , 

2kru sinc!>, 



= 0. 



(33) 



From the last line, we can see that ao remains constant 
as long as the contacts are non-sliding. This explains 
why ao = for all the MD simulations in Fig. ao = 
at the beginning of the simulation, and thus remains so 
until the disk rolls. This also explains why ao is constant 
when g is decreases in Figs. |S] and |H| 

Next let us calculate the disk's motion. Eqs. I|33[) can 
be solved for the velocities and integrated. In this way, 
the position r of the disk can be shown to be 





-mg 
r/(r 2 sin 2 . 



(34) 



Because no contact changes its status, the motion is per- 
fectly reversible. This is another way to see that the 
packing has no memory when tan</> < /i. 



2. Intermediate slopes: fi < ta,n<f> < 1/fi 

When tan^ > fi and t < tA, the contacts are sliding 
with T a = —fiR a and Tp — fiRp: 



( tan 2 



/caSc T = 2fccos 2 



V 



tan ( 



— ii tan i 


— fx tan 2 





Changes in contact status will occur, so it is necessary 
to subdivide the experiment into time intervals. The first 
change of status occurs at t = tA, when the derivative of 
the gravity changes sign, so we define our first interval 
to end at t = tA- We can find the particle displacements 
using the same method as before, and find, 



r y (t A ) 

As before, r x (t A ) = r e (i, 
a mfim , with 



mg n 



scc- 



2k 1 + fi tan cj) 
) = 0. We also have a (tA 



(36) 



tanc 



mg n 



fi 



1 + /i tan i 



(37) 



When the gravitational force begins to decrease, all 
contacts are non-sliding. The displacements must be cal- 
culated using the matrix given in Eq. (|32[1 . As before, 



the particle rises by a distance m(<? max — g*)/(2fc), while 
ao remains unchanged. The particle's position is now 



2k 



(38) 



It is important to realize that this is not the same as 
would have been obtained if the gravity was simply in- 
creased directly to g* . This is another expression of the 
packing's memory, but it is much more convenient to 
think of the memory being located in ao. Eq. I|38|) min- 
gles information about the past with information about 
the present, whereas Eq. (|37|l contains information only 
about the past. This is a consequence of using the ba- 
sis in Eq. JSJ, where what is determined directly by the 
imposed forces is separated from what is not. 

When the torque is applied, the contacts are still non- 
sliding, and they remain so until one of the contacts be- 
comes sliding. We can determine which contact slides 
first by checking if a mcm is greater than or less than a 2 . 
If 8mcm < 02, then the system meets the lower side of 
the funnel first, meaning contact a slides. Otherwise, it 
meets the upper edge first, and contact (3 slides. Using 
Eqs. I|17|) and (|37|l . we obtain that a mcm > &2 is equiva- 
lent to 



>g* 1 



fi 1 + tan 2 



tan < 



1 



(39) 



If this condition holds, then as r increases, ao decreases 
to maintain equality in Eq. (|9c|) . Eventually, equality is 
obtained in Eq. (|9a() . i.e., contact a begins to slip as well, 
and the disk rotates. This occurs when ao = 02, with 02 
given in Eq. 1)170. On the other hand, if Eq. (|39|l is not 
obeyed, than contact a slips first. Then ao increases 
until ao = a.2, and the disk begins to rotate. Evaluating 
Eq. (|3*9"|l for li = 0.5, <j) = 50°, one obtains <7 max ~ l-8g*, 
consistent with Fig. where simulations with g m & x /g = 
2, 4, 6, 8 all meet the upper edge of the funnel. 



3. Steep slopes: tan0 > 

For tan^ > l//i, all the calculation in the previous 
section applies. The only difference is that the slope of 
the line given by equality in Eq. (|9c|l changes sign. This 
line is shown in Fig. |SJ and if the system were to behave 
as before, it would follow this line upward indefinitely, 
leading to an infinite torque. However, the simulations 
show otherwise. When Eq. l|9"c|l is attained, the system 
suddenly jumps to the rotating state, in a behavior qual- 
itatively different from the other cases. This jump occurs 
when = —iiRp, so this suggests that we ought to cal- 
culate the displacements of the disk in this situation. We 
calculate the matrix A:aSc T as before, and assume that 
the torque is increasing. Calculating the velocity of the 
disk v in the usual way, one obtains 



2k (1 — fx tan 1 



[i — tan <j) + sec cf> esc < 

tan <j>{fi + tan 4>) I . (40) 
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Note the factor 1 — /ztan</> in the denominator, which is 
positive when tan cf> < l/fi and negative when tancj) > 
l//i. Let us pause for a moment to consider what this 
change of sign means for lu. Applying a positive torque 
to the disk means trying to rotate it in the positive di- 
rection, i.e., trying to impose u> > 0. When tan</> < l//i, 
Eq. H40fl shows that rotating the particle in this direc- 
tion generates contact forces that balance an increasing 
torque. Therefore, the configuration is stable: when r 
increases, the particle rotates in the positive direction, 
and the contact forces balance the increased torque. On 
the other hand, when tanc/> > l//i, the particle must be 
rotated in the negative direction to balance an increasing 
torque. But such a torque will always rotate the disk in 
the positive direction. Thus, the configuration is unsta- 
ble. The small increase in torque causes the particle to 
rotate in the positive direction as before. But this rota- 
tion causes the forces opposing the torque to decrease, 
which in turn causes the particle to rotate farther in the 
positive direction. This leads to a rapid change in the 
contact forces as the particle begins to rotate. Another 
way to see this is to calculate v T Qv using the value of v 
given above. One obtains 



v T Qv 



> sec' 



2kr 1 — /itan0 



(41) 



which becomes negative for tan<^> > l//x. Thus the par- 
ticle interactions no longer oppose the applied force, but 
amplify it. This leads to an exponential growth of the 
motion, and the rapid change in the forces observed in 
the simulations. 

We can now predict the yielding torque in the MD 
simulations. If Eq. 1)39(1 is not satisfied, the system meets 
the lower edge of the funnel that functions normally, and 
the disk rotates when r = T2 . An example of this is the 
MD simulation at <? max /5 ~ 2 in Fig. On the other 
hand, if Eq. I|39|) is satisfied, the system meets the other 
branch of the funnel, and immediately starts to rotate. 
We can calculate the torque where this happens because 
the value of ao is known. Thus inserting the value of ao 
given in Eq. (|37fl into Eq. (|9"c)l and solving for the torque 
gives 



r = mr sin ( 



9* +: 



tan <j) — fi fi tan 0—1 
tan 4> + fj, fj, tan tf> + 1 



This result accurately predicts the yielding torque, as 
shown in Fig. 



assumed to be rigid, the forces cannot be uniquely deter- 
mined. A certain combination of forces, whose amplitude 
is given by the coefficient ao in Eq. ©, causes no accel- 
eration, and thus cannot be deduced by examining the 
equations of static equilibrium. However, its value can 
be constrained by enforcing Coulomb friction and the ab- 
sence of cohesion at all the contacts. On the other hand, 
when the particles are assumed to be deformable, ao con- 
tains information about the history of the packing. This 
finding suggests that methods that average over different 
possible force networks Q are really averaging over past 
histories of the packing. 

Secondly, we were able to precisely understand the con- 
nection between indeterminacy and motion. Depending 
on the angle </>, the onset of motion of the disk can take 
two different forms. When tan</> < l/fi, the range of al- 
lowed values of ao decreases continuously, and vanishes 
precisely at the onset of motion, and the torque needed 
to turn the disk is independent of method and history 
of the packing. Furthermore, the forces change slowly 
and continuously. A second mode of failure is more dra- 
matic. When tan<^ > l/fi, the indeterminacy is infinite, 
i.e., there is no upper bound on ao. At the onset of mo- 
tion, the forces change discontinuously (on the time scale 
of changes in the imposed torque) , and the torque needed 
to turn the disk depends on the history and method. This 
mode of failure is related to a curious situation in which 
the contact forces no longer oppose the applied forces, 
but enhance it. 

This work shows that giving the CD iterative solver 
the the solution of the last time step as the initial guess 
is not just a numerical trick to accelerate convergence. It 
is an essential part of the method that represents a real 
physical process at work in granular material, namely the 
memory of the packing encoded in the combinations of 
force that cause no acceleration. In certain situations, 
this memory could have an important effect on the be- 
havior of the packing. 

It is hoped that the findings presented in this paper 
will deepen and clarify our understanding of the quasi- 
static deformation of much larger packings. When the 
number of particles is larger, there is not just one unde- 
termined combination of forces, but many. However, one 
can expect that the effects studied here will be present 
in these more complicated situations, perhaps combined 
with each other in new and unexpected ways. 
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